One of the possible methods of stabilizing the Stokes equations is the DG (Discontinuous Galerkin) method. In this method, the basis functions used to approximate and test the solutions are "cut" into pieces with scissors at the boundary of individual elements. So if we want to use B-spline polynomials, we can continue to use them, but over each element they will be cut into local pieces.
In this chapter, we build on the DG formulation of Stokes equations [1].
We use the formulation derived in the previous chapter.
Additionally, the average of the values of the dissected basis functions on the edges of two adjacent elements is introduced
\( a({\bf u},{\bf v}) + b(p,{\bf v }) = f({\bf v}) \quad \forall {\bf v} \)
\( b(q,{\bf u } ) = 0 \quad \forall q \)
\( a({\bf u}, {\bf v}) =\int_{\Omega} \nabla {\bf u} : \nabla {\bf v }dxdy = \int_{\Omega} \frac{\partial u_1}{\partial x}\frac{\partial v_1}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_1}{\partial y}\frac{\partial v_1}{\partial y}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial x}\frac{\partial v_2}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial y }\frac{\partial v_2}{\partial y }dxdy \)
\( b(p, {\bf v}) = \int_{\Omega} p div {\bf v} dxdy= \int_{\Omega}p \frac{\partial v_1}{\partial x} dxdy+ \int_{\Omega} p \frac{\partial v_2}{\partial y } dxdy \)
\( f({\bf v} )= \int_{\Omega} \frac{\partial h(x,y)}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial h(x,y)}{\partial y}\frac{\partial v_2(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial h(x,y)}{\partial n } v_1(x,y) dS \)
In the case of formulation for the DG method, it is necessary to introduce additional terms describing jumps of basis functions \( [\psi(x)]=\psi|_{T_1}(x)-\psi|_{T_2}(x) \) where \( T_1 \) i \( T_2 \) stands for two adjacent elements, and we calculate here the difference from the "cleaved" value of the B-spline function marked \( \psi \) on the element on one side and on the other side.
Additionally, the average of the values of the dissected base functions on the edges of two adjacent elements is entered \( \{\psi\} = \frac{\psi|_{T_1}(x)+\psi|_{T_2}(x)}{2} \).
Our functional \( a({\bf u },{\bf v }) \) should be modified for DG method as follows
\( a_{DG}({\bf u }, {\bf v }) = \sum_{K \in \Omega_h } \int_K \left( \frac{\partial u_1 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_1 }{\partial y } \frac{\partial v_1 }{\partial y } + \frac{\partial u_2 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_2 }{\partial y } \frac{\partial v_1 }{\partial y } \right) \)
\( + \sum_{F\in F_h} \int_F \frac{\eta}{h_F} \left( [\frac{\partial u_1}{\partial x} ] [\frac{\partial v_1}{\partial x}] + [\frac{\partial u_1}{\partial y}] [\frac{\partial v_1}{\partial y }] + [\frac{\partial u_2}{\partial x} ] [\frac{\partial v_2}{\partial x}] + [\frac{\partial u_2}{\partial y}] [\frac{\partial v_2}{\partial y }]\right) \)
\( -\sum_{F\in F_h} \int_F \left( \{\frac{\partial u_1}{\partial n_1} \} [\frac{\partial v_1}{\partial x}] + \{ \frac{\partial u_1}{\partial n_2} \} [\frac{\partial v_1}{\partial y } ] + \{\frac{\partial u_2}{\partial n_1} \} [\frac{\partial v_2}{\partial x}] + \{ \frac{\partial u_2}{\partial n_2} \} [\frac{\partial v_2}{\partial y } ]\right) \)
\( - \sum_{F\in F_h} \int_F \left( [\frac{\partial u_1}{\partial x} ] \{ \frac{\partial v_1}{\partial n_1} \} + [\frac{\partial u_1}{\partial y} ] \{\frac{\partial v_1}{\partial n_2 } \} + [\frac{\partial u_2}{\partial x} ] \{ \frac{\partial v_2}{\partial n_1} \} + [\frac{\partial u_2}{\partial y} ] \{\frac{\partial v_2}{\partial n_2 } \} \right) \)
In the above formulas, the integrals are computed element by element. Thus \( \Omega_h \) means the collection of all elements from the patch on which our B-spline functions are spread, \( K \in \Omega_h \)denotes a single element selected from a family of elements (so we sum our element integrals over all elements in the patch), \( F_h \) stands for the collection of all edges between patch elements, \( F\in F_h \) denotes a single edge selected from this family (so we sum our edge integrals over all edges from the patch), \( h_F \) denotes the edge length \( F \), and \( \eta \) denotes the stabilization parameter. Its proper value needs to be found experimentally.
Unfortunately, the disadvantage of this stabilization method is the fact that we have to run a few calculations and find this parameter value \( \eta \) for which the simulation will give correct results.Typically, you experiment with values that differ by an order of magnitude, e.g. 0.1, 1, 10, 100, 1000, 10000 etc.
Symbol \( {\bf n }=(n_1,n_2) \) denotes the unit vector perpendicular to the edge \( F \). In our case, when we consider the B-spline base functions spanned over a regular 2D patch of elements, we assume \( {\bf n }=(1,0) \) or \( {\bf n }=(0,1) \).
Our functional \( b(p ,{\bf v }) \) should be modified for DG method as follows
\( b_{DG}(p, {\bf v})= - \sum_{K\in \Omega_h } \int_K \left ( p \frac{\partial v_1}{\partial x}+ p \frac{\partial v_2}{\partial y } \right) \)
\( + \sum_{F \in F_h } \int_F\left( [v_1] n_1 p+ [v_2] n_2 p \right) \)
The meaning of the symbols in integrals is similar to that of a functional \( a({\bf u },{\bf v }) \).
Our stabilized Stokes problem now looks like this:
\( a_{DG}({\bf u},{\bf v}) + b_{DG}(p,{\bf v }) = f({\bf v }) \quad \forall {\bf v} \)
\( b_{DG}(q,{\bf u } ) + {\color{red} {S_{DG}(p,q) }}= 0 \quad \forall q \)
According to the book Daniele Antoni Di Pietro, Alexandre Ern, Mathematical Aspects of Discontinuous Galerkin Methods, to stabilize the Stokes problem we need to add the term
\( {\color{red}{S_{DG}(p,q)}} = \sum_{F\in F_h^{internal} } h_F \int_F [p][q] \)
where \( F_h^{internal} \) stands for all inner edges, so we sum here the integrals calculated on all edges of the elements inside the patch (avoid the edges located on the edge). As before, \( h_F \) denotes the edge length.
In order to solve the Stokes equations stabilized by the DG method, we generate a similar system of equations, taking into account the above modifications.
It is also possible to consider various modifications.
We will assume the following facts about the approximating and testing functions. Let us assume that the test functions will be continuous, in the sense that over individual segments of the knot vector in the isogeometric analysis, or over individual finite elements in the sense of the classical finite element method, they will be continuous, but will be discontinuous between the elements or between the knot vector segments. So we have "broken" test polynomials. Let us also assume that approximation functions are continuous. Then, all terms containing jumps and averages from the approximating function will disappear, in particular we will get modified expressions
\( a_{DG}({\bf u }, {\bf v }) = \sum_{K \in \Omega_h } \int_K \left( \frac{\partial u_1 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_1 }{\partial y } \frac{\partial v_1 }{\partial y } + \frac{\partial u_2 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_2 }{\partial y } \frac{\partial v_1 }{\partial y } \right) \)
\( -\sum_{F\in F_h} \int_F \left( \{\frac{\partial u_1}{\partial n_1} \} [\frac{\partial v_1}{\partial x}] + \{ \frac{\partial u_1}{\partial n_2} \} [\frac{\partial v_1}{\partial y } ] + \{\frac{\partial u_2}{\partial n_1} \} [\frac{\partial v_2}{\partial x}] + \{ \frac{\partial u_2}{\partial n_2} \} [\frac{\partial v_2}{\partial y } ]\right) \)
\( b_{DG}(p, {\bf v})= - \sum_{K\in \Omega_h } \int_K \left ( p \frac{\partial v_1}{\partial x}+ p \frac{\partial v_2}{\partial y } \right) \)
\( + \sum_{F \in F_h } \int_F\left( [v_1] n_1 p+ [v_2] n_2 p \right) \)
\( {\color{red}{S_{DG}(p,q)}} = 0 \)
We will assume the following facts about the approximating and testing functions. Let us assume that the approximating functions will be continuous pieces in the sense that over individual segments of the knot vector in the isogeometric analysis, or over individual finite elements in the sense of the classical finite element method, they will be continuous, but will be discontinuous between the elements or between the segments of the knot vector. So we have "broken" approximation polynomials. Let us also assume that the test functions are continuous. Then we have
\( a_{DG}({\bf u }, {\bf v }) = \sum_{K \in \Omega_h } \int_K \left( \frac{\partial u_1 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_1 }{\partial y } \frac{\partial v_1 }{\partial y } + \frac{\partial u_2 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_2 }{\partial y } \frac{\partial v_1 }{\partial y } \right) \)
\( - \sum_{F\in F_h} \int_F \left( [\frac{\partial u_1}{\partial x} ] \{ \frac{\partial v_1}{\partial n_1} \} + [\frac{\partial u_1}{\partial y} ] \{\frac{\partial v_1}{\partial n_2 } \} + [\frac{\partial u_2}{\partial x} ] \{ \frac{\partial v_2}{\partial n_1} \} + [\frac{\partial u_2}{\partial y} ] \{\frac{\partial v_2}{\partial n_2 } \} \right) \)
\( b_{DG}(p, {\bf v})= - \sum_{K\in \Omega_h } \int_K \left ( p \frac{\partial v_1}{\partial x}+ p \frac{\partial v_2}{\partial y } \right) \)
\( {\color{red}{S_{DG}(p,q)}} = 0 \)
If we put continuous approximation and test spaces into the DG formulation, we will get the same original Stokes formulation without stabilization. Different variants of approximation and testing spaces can be helpful when testing various forms of stabilization and for debugging the code generating the matrix of the DG method.